Fast simulations of charges moving in response to varying dielectric media 
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For charged systems in heterogeneous dielectric media, a key obstacle for fast molecular dynamics (MD) 
simulations is the need to solve the Poisson equation in the media. This obstacle can be bypassed using MD 
methods that treat the local polarization charge density as a dynamic variable, but such approaches require 
access to a true free energy functional; which evaluates to the equilibrium electrostatic energy at its minimum. 
In this letter, we derive the needed functional. As an application, we develop a Car-Parrinello MD method for 
the simulation of free charges present near a spherical emulsion droplet separating two immiscible liquids with 
different dielectric constants. Our results show the presence of non-monotonic ionic profiles in the dielectric 
with lower dielectric constant. 



It is hard to overstate the importance of the electrostatic 
force in biological and soft-matter sciences. Electrostatic in- 
teractions play a major role in determining the structure and 
function of several biological macromolecules, such as pro- 
teins and DNA 01 12]. In cell signaling, the creation of elec- 
trical potential differences and transfer of ions are of chief 
importance . On the other hand, electrostatic forces allow 
the stabilization of many synthetic structures, endowed with 
remarkable properties: self-assembled colloidal dispersions 
01, overcharged surfaces 0], patterned surfaces by competi- 
tion between short-range and Coulombic interactions @], and 
faceted thin shells 0], to name a few. 

Investigation of biological and soft-matter systems there- 
fore requires an accurate consideration of electrostatic inter- 
actions. In many situations, computational methods that treat 
ions individually are necessary. This is the case, for exam- 
ple, where finite size effects become significant, where the 
medium has a complex geometry, or where the dielectric re- 
sponse of the medium is inhomogeneous. Computing prop- 
erties of such systems via numerical simulation involves its 
own challenges. Due to the long range of the Coulomb force, 
a system of N charges requires an expensive 0(N 2 ) force 
calculation at every simulation step. Attempts to ameliorate 
this scaling have resulted in the development of several meth- 
ods: Ewald summation, particle-mesh methods, fast m ultip ole 
methods jfj|, and the local electrostatics algorithms 1 ill . In 
this Letter, we focus on the problems associated with the pres- 
ence of dielectric heterogeneities in the medium. 

The presence of free charges in a medium polarizes its un- 
charged constituents, leading to a complex behavior for the 
electric field and the polarization field itself. Accurate inves- 
tigations of systems with electrostatic interactions should in- 
corporate this dielectric response of the medium. In many 
cases, it is sufficient to consider the polarization effects by de- 
scribing the dielectric properties with a spatially varying di- 
electric constant. In the simplest case of a uniform dielectric 
response, a single dielectric constant enters the coarse grained 
model, and the simulation proceeds as it would in free space, 
albeit with a scaled Coulomb's law. However, most real sit- 
uations involve regions with different dielectric response, as 
is the case for proteins within an aqueous cellular medium or 



for emulsions where oil and water are partitioned fl2)l ■ In the 
presence of this varying dielectric response, the simplest form 
of Coulomb's law breaks down and one has to solve the Pois- 
son equation, at each simulation step, to obtain the necessary 
force information for the propagation of charges. This ad- 
versely affects the stability and efficiency of the resulting nu- 
merical procedure. Because of these challenges, the problem 
of treating variable dielectric response in c harg e simulations 
continues to receive much attention ifl^ [T3U2OII . 

In this Letter, we present a variational formulation of elec- 
trostatics that is applicable to problems involving dielectric 
heterogeneities. We construct an energy functional with the 
polarization charge density as the sole variational field. This 
functional works for arbitrary free charge configurations and 
any kind of spatial variation in dielectric response. As we 
review later, these characteristics have not been realized in 
previously proposed functionals. We explicitly specialize the 
functional to the case of sharp interfaces, where only the 
surface polarization charge density needs to be extremized. 
We also demonstrate that our functional can be used in sim- 
ulations of charged systems by employing a Car-Parrniello 
molecular dynamics (CPMD) scheme 12 ill where the surface 
polarization charge density is treated as a dynamical variable. 

The polarization charge density uj is defined by the rela- 
tion oj = — V • P, where P is the polarization field. We 
assume that the medium polarization obeys linear response, 
P = \Ei, where \ is the susceptibility and E is the electric 
field. Employing the notation p for the free charge density 



and G(r,r') = |r 
functional reads: 
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for the bare Green's function, our 
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where Q r [uj] is both a functional of w(r) and a function of r, 
and is defined as 
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Extremization of ^\oj] leads to the equality: u)(r) = fi r 
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which is the correct physical relation that lo must satisfy. In 
spite of the complex dependence on w(r), our functional re- 
tains a simple interpretation at equilibrium as, owing to the ex- 
tremum condition, its second term vanishes and the first term 
becomes the true electrostatic energy U = \ J p(r)(/)(r)d 3 r, 
where <p(r) is the electrostatic potential. Furthermore, it can 
be shown that the second variation of ,^[lo] is positive, imply- 
ing ,^[lo] becomes a minimum at its extremum. We give the 
proofs elsewhere 12211 . 

A variety of functionals employing various electrostatic 
quantities as field variables have been proposed Ull [l4l ITtI 
2^1. Many of these functionals are not energy functionals 



23 



riy o l 



231 126 , 127H : n amely, they either become a maximum at 



extremum 1 26ll27l| or evaluate to negative electrostatic energy 
at equilibrium Il4, 12311 . This rules out their use in dynamical 
optimization methods 1 21 ] . Attard lfl7ll has given an energy 
functional of the surface polarization charge density, but his 
functional is derived for a specific system that involves all free 
charges to be constrained in one uniform dielectric medium. 
Other energy functionals ll U l24l 12511 involve relatively ex- 
pensive vector function variables, requiring three-dimensional 
vectorial specification 1 1 il l 1 311 . In problems where the di- 



electric response can be modeled as piecewise uniform, our 
functional reduces to a functional of the surface polarization 
charge density, which requires only a two-dimensional scalar 
specification and offers distinct numerical advantages over the 
vector variables. 

Our main result, Eq. (Q]i, can be derived as follows. We 
begin with the standard expression for the electrostatic en- 
ergy written as a functional: & [E] = J e (r) |E ( r)\ d 3 r, 
where e is the dielectric permittivity. Next, following 12911 . we 
include Gauss's law as a constraint to this functional via the 
Lagrange multiplier <fi: 
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Note that we take & to depend parametrically on p. Using 
e = 1 + 47rx and P = xE, we introduce P in (O to obtain 
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Variations of with respect to E and <f> allow us to eliminate 
all variables in favor of P, leading to 
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[P] has been derived before using different derivations 
|24j,|25|]. It can be shown that & [P] is an energy functional; 
that is, its minimum computes the equilibrium electrostatic 
energy lfl7h . We wish to transform J? [P] to an energy func- 
tional of lo. This transition begins by inserting the definition 



of a; in © via a Lagrange multiplier tp: 

^[P,W,V] = J l ^f d3r + l J J (Pr + Wr)Gr,r< 

x (p r , +uj r >) d 3 r'd 3 r - J ip T (lo t + V • P r ) d 3 r. 
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We note that ip can be shown to coincide with the electrostatic 
potential <fi at equilibrium. Taking variations of the above 
functional with respect to lo, P, and ijj gives: 

V>(r) = J G T y{p T ,+u r ,)d 3 r', (7) 

P r = -XrV r J G r . r , ( Pr , + Ur >) d 3 r' , (8) 

uj r = V • (xvV J G r y (j)^ + uj r ,)d 3 r'^j . (9) 

It is evident from (0 and ((H) that ip is the electrostatic poten- 
tial. In addition, following (O, we can recast Eq. (O as the 
equality of two quantities: the polarization charge density lo 
and the operator O [lo] . This equality is precisely the extremum 
condition for cf\io\ Substitution of ip, P, and lo from Q, (|S), 
and (O respectively, into (|6]l leads to our central result, the 
functional in (Q]). 

We note that not all substitutions lead to our result. 
For example, eliminating P and tp from © using (0 and 
(© gives a functional I[lo] with the functional density: 
PtGt,t' (pr' + ^r'[ w ]) /2 — w r G T y (to T ' — Q r '[uj]) /2. Upon 
extremization, I[lo] singles out the correct physical quantity, 
but becomes a maximum at equilibrium. In fact, I[lo] is ex- 
actly the negative of the functional in 11411 . neither of which 
are energy functionals. Functionals ,^[lo] and I[lo] share a 
common structure: the expression for the total electrostatic 
energy (the first term in either functionals) is constrained 
by the correct physical relation that lo must satisfy, namely 
lo — Qj[lo] =0. The only but crucial difference between these 
two functionals is the choice of the Lagrange multiplier that 
enforces this constraint. In past attempts, the Lagrange multi- 
plier that leads to an energy functional of lo remained elusive. 
Our current formulation finds the desired multiplier. 

We now consider the application of our functional to the 
problem of point charges in a system with piecewise-uniform 
dielectrics. For the sake of brevity, we will restrict ourselves 
to two uniform dielectrics, with different permittivities ei and 
£2, separated by a single sharp interface. Extension to mul- 
tiple dielectrics and interfaces is straightforward. We assume 
that free charges reside only in the bulk of the dielectric. The 

Ejli <H S ( r - r *)' where 



(5) free charge density is p(r) 



is the charge of the i particle and N is the total number of 
charges. For this system, the induced charge density in the 
bulk has the form Wbuik(r) = (l/ e ( r ) — 1) P ( r )' which leads 
to an effective charge density of p(r)/e(r). Also, the gradient 
of x vanishes everywhere except at the interface. Therefore 
several volume integrals in Eq. (fTJ reduce to surface integrals 
and the original functional ,jF[lo] is transformed to a functional 
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of the surface induced charge density: 
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where uj(s) is the induced charge density at the position s on 

the interface, and K , K , and K are, respectively, the effective 
potentials of interaction between two free charges, between a 
free charge and an induced charge, and between two induced 
charges. These effective interactions are given by: 

1 1_1 1=1 

K r r i — — Gr,r' ~\~ — Gr,r' ~F" — Gr.r' 



Gs,v — (2e m — 1) G r ,s 2G r 



K s , s , = e rn (e m - 1) G SjS , - (2e m - 1) G s , s , + G^,. (11) 

Here e TO = (ei + £2)/2 is the permittivity at the interface and 
functions G and G are defined as: 



G a ,b = £d / G a , u n u ■ VG u ,b rf 2 " 

G a ,b = e 2 d n u - VG a , u G u . v n v ■ VG v ,b d 2 ud 2 v, (12) 



where = |e2 — ei|/47r is the permittivity jump at the in- 
terface, a, b are arbitrary position vectors, u, v are position 
vectors of arbitrary interfacial points, and h is the unit normal 
to the interface taken to point in the direction of increasing 
permittivity. 

In addition to providing a complete reformulation of elec- 
trostatics in heterogeneous media, our formalism has immedi- 
ate applications to important practical problems. Since 
is an energy functional, it can be used for simulating free 
charges in heterogeneous media which, as described above, 
are basic models for phenomena in both biological and syn- 
thetic settings. The simplest simulation schemes 00 for 
these systems require, in some way, the solution of the ex- 
tremum condition, u — = 0, at each step. However, 

when an energy functional is available, new approaches are 
possible, such as the use of CPMD method 12 ill . In this ap- 
proach, u> is treated as a dynamical variable and is assigned 
a mass. The dynamic equations for the system follow from a 
Lagrangian that contains an additional kinetic energy term for 
u). The kinetic term is constructed so that lu remains close to 
the exact polarization charge distribution at all times. In other 
words, we replace the expensive solution of the Poisson equa- 
tion at each simulation step with an on-the-fly computation of 
the polarization effects. 

We apply our method to free charges in piecewise-uniform 
dielectrics, where the interfaces between the different uniform 
dielectrics are closed surfaces of finite area. For simplicity, we 
only consider impenetrable boundaries such that each region 
has a fixed set of ions. We partition the surface boundaries 



into M finite elements, and to each element k we assign an 
average induced charge density LUk and a fictitious mass fik- 
For the system with N free ions, we can write a Lagrangian 
for the extended system of TV + M particles as: 



N 1 M 1 



J^r,] - Jf[r,]. (13) 



fc=i 



The first term is the kinetic energy of N ions with masses rrii 
and position r;. The second is a fictitious kinetic energy for 
the surface induced charge density. The electrostatic potential 
energy of the system computed by using our functional con- 
stitutes the third term. And the final term contains a set of 
truncated Lennard-Jones potentials to model the hard-core of 
the particles and the impenetrability of the surfaces. 

Starting from a point in the extended configuration space, 
we generate its dynamics via standard MD technique, using 
the equations of motion derived from the Lagrangian £ for 
the ions and the fictitious induced charge values. To simu- 
late the behavior of the physical system at finite temperature 
T, we couple the augmented system to a set of Nose-Hoover 
thermostats. The ions couple to a thermostat at temperature 
T, while the induced charge values couple to one at much 
lower temperature T2. This allows the evaluated energy of the 
physical system to be close to its thermal equilibrium value by 
limiting the contribution of the fictitious kinetic energy. This 
two-temperature approach is a standard feature of CPMD ll30l - 
3211 . The masses of the induced charge degrees of freedom are 
chosen so as to make their energy contribution small. In prac- 
tice, we choose these masses to be proportional to the areas of 
the finite elements, and the proportionality constant is chosen 
to optimize the stability of the simulation. Another feature of 
the system we simulate is that, as a result of Gauss's law, the 
net induced charge at each interface is a constant. In our sim- 
ulations, we choose to directly enforce this constraint at each 
step via the SHAKE-RATTLE routine J33]. 

We have used the CPMD method outlined above to simu- 
late ions near a spherical interface separating media of differ- 
ent permittivities. The simulations have been carried out for 
various values of permittivities, ion concentrations and ion va- 
lencies. As a test case, we considered a model for charged col- 
loidal dispersions, where mobile charges are present in only 
one of the two dielectrics. We obtained results that match 
those previously published \ \5, 2(J. In this Letter, we focus 



on a much less explored problem of ions present in both sides 
of the spherical interface. This system can be considered as a 
model for a li quid -liquid emulsion droplet in the presence of 
an electrolyte 11341. 13511 . We consider the case where the ions 
do not cross the interface, just like in the experiments of l36ll . 
We model the ions as repulsive Lennard-Jones (LJ) spheres 
of diameter a = 0.357 nm. The spherical interface of radius 
a = lOcr separates the two media: the interior dielectric has 
permittivity e m = 80, while the exterior dielectric has per- 
mittivity e out = 35. The whole system of ions and dielectric 
media is taken to be in a spherical simulation cell of diameter 
b = 20a such that the centers of the two spheres coincide, 
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tion is monotonic and gets more pronounced for higher values 
of c; n . On the other hand, in the exterior, lower-permittivity 
dielectric, ions prefer to accumulate near the interface (see in- 
set in Fig. [TJ. We also see that ionic profiles on this side of 
the interface are non-monotonic. This is due to a combination 



FIG. 1 . Ionic density profiles for different concentrations of monova- 
lent electrolyte inside the dielectric sphere. The outside salt concen- 
tration is held fixed at 0.3 M for all curves. The inset shows the ac- 
cumulation of ions near the interface on the lower permittivity side. 
The dashed line in the middle shows the position of the interface. 
On the left side we sketch our model for the liquid-liquid emulsion 
droplet in the presence of a monovalent electrolyte (not to scale). 



see the sketch in Fig. [I] Both the interface and the simulation 
cell boundary are modeled as repulsive LJ walls. We consider 
monovalent electrolyte (1:1) at T = 298 K and c m (c out ) de- 
note the salt concentrations inside (outside) the spherical in- 
terface. The interface is discretized with nearly M = 2000 
points, and the parameters associated with the CPMD simula- 
tion are: A = 0.001, \i k ~ 5 - 10, T 2 = 0.001T. Here, A is 
the simulation timestep and these values are given in LJ units 
(energy: fc^T, length: a). 

Our simulation results pass two tests of stability and accu- 
racy. We have analyzed the energy of the simulated system as 
a whole as well as the energies of its subsystems. Fluctuations 
in the total energy of the physical system were found to be 50 
times smaller than those in the physical kinetic energy, imply- 
ing good energy conservation. Also, the energy of the ficti- 
tious kinetic modes was kept very close to zero at all times. 
Our second test relates to the effectiveness of our scheme to 
reproduce the correct polarization charge distribution. At reg- 
ular intervals during the course of the simulation of a number 
of specific cases, we collected and stored the ion coordinates 
and surface charge densities. Then, we carried out an ordinary 
minimization of the functional to determine the exact induced 
density, and compared it to the distributions obtained in the 
simulation. Our on-the-fly method results were within 2% of 
those obtained with direct minimization. 

Our simulations of the 1:1 electrolyte lead to the density 
profiles shown in Fig. 03 We consider different values of Ci n 
while maintaining c out at 0.3 M, thus conducting a study sim- 
ilar in spirit to the experiment in lt36ll . The ion distributions, 
for all concentrations, reach a constant value in the bulk on 
either side but show interesting features near the interface. On 
the side having the higher value of permittivity, the ion den- 
sity is depleted near the interface, which is largely the result 
of the repulsion due to the induced surface charge. The deple- 



of Coulombic depletion near hard wall B7U38I1 and attractive 
surface polarization charge effects. We further observe that 
increasing the internal salt concentration enhances the accu- 
mulation of external ions near the interface. Several of these 
features have been previously observed for planar interfaces 
lfl6ll and attributed generally to the same basic reasons. 

We have presented the solution to the long-standing prob- 
lem of providing a true free energy functional for electrostat- 
ics that employs the polarization charge density as the vari- 
ational field. This formulation has many applications, and 
we have used it to develop an efficient CPMD simulation for 
a set of point charges present in two dielectric media. Fur- 
ther applications of this approach, and its coupling with other 
techniques such as Ewald summation, will produce consider- 
able gains in computational efficiency over currently available 
methods. As the functional is general for linear media, simu- 
lation approaches derived from it can be applied to many other 
systems, such as those with arbitrary interface shapes or mov- 
ing boundaries. 
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